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Abstract. In this contribution, we introduce numerical continuation methods and 
bifurcation theory, techniques which find their roots in the study of dynamical systems, 
to the problem of tracing the parameter dependence of bound and resonant states of the 
quantum mechanical Schrodinger equation. We extend previous work on the subject [ ] 
to systems of coupled equations. 

Bound and resonant states of the Schrodinger equation can be determined through the 
poles of the S-matrix, a quantity that can be derived from the asymptotic form of the 
wave fimction. We introduce a regularization procedure that essentially transforms 
the S-matrix into its inverse and improves its smoothness properties, thus making it 
amenable to numerical continuation. This allows us to automate the process of tracking 
boimd and resonant states when parameters in the Schrodinger equation are varied. 
We have applied this approach to a number of model problems with satisfying results. 
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1 Introduction 

The appearance of resonances is of ever-growing interest in the study of wave phenomena 
as they are considered among the most important features of systems described by wave 
equations. They appear in systems that are penetrable by an impacting wave. Such systems 
allow the interior field to couple to the external domain which leaves a characteristic 
fingerprint on the far-field pattern of the scattered wave. Many examples appear naturally 
in acoustic scattering [ ] and in fluid-mechanical structure interaction [ ]. In all cases the 
appearance of resonant states has a profound and important influence on the system's 
dynamics. 

In the context of quantum mechanical scattering, resonant behavior also strongly 
influences the interactions between microscopic particles, which in turn has its influence 
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on the reactivity of molecules and atoms described by such quantum mechanical models. 
In molecular systems these resonances can easily turn into bound states if the molecular 
configuration changes. 

In [1] we developed a framework for applying numerical continuation techniques in 
the context of bound states and resonances in spherically symmetric short-range potentials. 
We have shown that numerical continuation methods, originally developed in the study 
of dynamical systems, can be applied successfully to track bound and resonant states 
efficiently in terms of a varying system parameter. Moreover, this technique can be used 
to reveal subtle and interesting transitions and connections between states automatically. 

The present work focuses on the extension of that procedure to coupled channel short- 
range systems. This extension is a logical step towards automated, efficient and robust 
methods for the study of interactions in scattering experiments. In all generality, these 
techniques can be applied to systems of coupled Helmholtz equations with variable wave 
numbers, as long as the short-range conditions are met. 

The outline of the paper is as follows. Section 2 sets the coupled channel Schrodinger 
equations in the context of non spherically symmetric quantum mechanical problems. 
In section 3 a regularization procedure is discussed that allows application of numerical 
continuation even though the underlying functions that characterize resonances and 
bound states are numerically and analytically not very well-suited. Section 4 provides a 
brief overview of basic numerical continuation methods and gives some pointers on the 
available implementations. Finally, in section 5 we present several results obtained with 
our implementation of the discussed methods. 

2 Quantum scattering in coupled channel problems 

In this paper we discuss a coupled channel problem that derives from a one particle 
Schrodinger equation with a non-spherical potential 



where A is the three-dimensional Laplacian, V{r,\) is a potential with a system parameter 
A and E is the complex-valued energy of the system. The problem is such that for 
all r outside a bounded domain F(r,A) ^ 0, i.e. the potential becomes negligibly small. 
Formally, the limitation to potentials that are negligible outside a certain radius is termed 
as the restriction to so called short-range potentials: V{r,\) must decay faster than r^^ as 
r = |r I — 7- 00 and must be less singular than r^'^ in the origin r = |r | — )■ [ ]. The long-range 
Coulomb interaction requires a significantly different approach and is not discussed here. 
The boundary conditions that are appropriate in the short-range case require the solution 
to be zero at the origin and force the solution to tend to a linear combination of free waves 
(i.e. solutions of (2.1) for y = 0) at infinity [ ]. 

In this section we briefly introduce elements from quantum scattering and partial wave 
expansion to arrive at the concept of resonant states. In the subsequent sections we then 
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focus on applying numerical continuation to study the dependence of such resonances on 
the system parameter A in (2.1). For this reason, we will faithfully record the A dependence 
in our notations, even if it is at times somewhat cumbersome. 

Equation (2.1) can be written in spherical coordinates {r,6,q)) around the center of 
the system. The differential operator A then splits into angular and radial differential 
operators and the solution can be expanded as a sum 

CO / 

rp{r) = xp{rAcp) = Y: L ^Ur)YT ,<?) , (2.2) 

/=0m=-/ 

where Yl"{6,q)) are the spherical harmonics, eigenf unctions of the angular differential 
operators in the Laplacian. In physics this decomposition is commonly referred to as 
the partial wave expansion. The labels {l,m) are intimately connected to the irreducible 
representations of the rotational symmetry groups SO(3) D SO(2) of equation (2.1). Each 
(/,m) -component of (2.2) is known as a partial wave and the system is said to be modeled 
by multiple partial wave channels. We are now interested in localizing the resonances in 
such systems. 



2.1 Single-channel scattering 

To introduce some notations and concepts we first briefly look at a spherically symmetric 
potential, i.e. V{r,\) = V{r,A). In this case, the partial wave channels are decoupled and 
the radial wave function xpi^ is identical for all m. Hence we can drop the index m and 
equation (2.1) turns into: 



1 /(/ + 1) 



- ipi (r, A) + V{r,A)ipi (r, A) = 0, (2.3) 



for each /. The boundary condition at r = specifies xpi{r,\) =0 and at large r the solution 
must be a linear combination of the spherical Riccati-Hankel functions (see appendix A), 
the free incoming and outgoing waves: 

rp,{r,A)^'-{h-{kr)-h+{kr)Si{K\)y (2.4) 

The spherical Riccati-Hankel functions are the solutions obtained in the absence of the 
potential term in (2.3), i.e. when the equation reduces to a Helmholtz equation with a 
constant wave number k = y^lfiE. The S/ (A:, A) is called the S-matrix for channel / and it 
determines the scattering properties associated with potential V{r,\). 



2.2 Multi-channel scattering 

For a non-spherical potential the partial wave channels do not decouple. The wave 
function must be represented by a sum as in equation (2.2), which is typically truncated 
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at an Zmax- In most cases of physical interest the potential still has axial symmetry. It is 
well known that as a consequence the channels with different m are decoupled and the 
solutions can be represented by 

^ max 

^P{r) = Y,%m{r)Yr{e,q>). (2.5) 

1=0 

Upon substitution in (2.1), this generates a separate set of equations for each m. For the 
purpose of our exposition we may, without loss of generality, take m = and drop the 
index m in the notation for the radial wave function 

^max 

ip{r) = Y2Mr)y?{0,<p). (2.6) 

1=0 

This significantly simplifies notations further on in the discussion. Lifting the restriction of 
axial symmetry presents no fundamental difficulty: the number of channels increases and 
summations over m have to be reintroduced. We will however consider only the axially 
symmetric case with m = as it provides a sufficient basis to demonstrate that numerical 
continuation can be used to track resonant states in multi-channel systems. 

Insertion of (2.6) in (2.1), followed by a projection onto each of the partial wave 
channels (by multiplying with Y^,, integrating over the angular variables and using the 
orthonormality of the spherical harmonics) leads to a set of coupled radial equations 

/ 1 j2 1(1-1-1) \ 'max 

(-2^*^ + ^-^) «"<''^' + LV„.(r,m-M = 0, (2.7) 
for each channel / and where 

y„,(r,A) = Jsm9d9d(pYf{e,(p)V{r,e,A)Yj',{e,(p), (2.8) 

and where * denotes the complex conjugate. For every partial wave tpi the boundary 
condition at r = is still j/j/(r,A) =0 while for large r one formulates /max + 1 independent 
boundary conditions (one for each /') 

xpl'{r,\) ^ (h^{kr)Sjr-hl{kr)Si,{k,\)), (2.9) 

where Sw is the usual Kronecker symbol. 

Each boundary condition corresponds to an incoming wave in a different channel /' 
and gives rise to an independent associated solution with its channel components labeled 
by /. A key difference with the spherical case is that, due to the non-spherical action of 
the potential, outgoing wave components appear in channels other than the incoming 
channel, as indicated by the S-matrix defined by its elements Sw (A:, A). 
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We now rewrite the above in terms of matrices by introducing the square matrices T, 
V and the identity matrix I with dimensions equal to the number of channels: 



1 1(1 + 1] 



(y(r,A))„, = y„(r,A). 



(2.10) 



Eq. (2.7) then is equivalent to 



(T(r)-£I+y(r,A))Y''(r,A)=0, 



(2.11) 



where the wave function is a vector function with dimension equal to the number of 
channels, corresponding to a solution with (2.9) (easily rewitten in vector format). In 
multi-channel calculations it is however customary to introduce the square matrix T 
whose columns are the f ' (corresponding to (2.9)) for each of the This then leads to 
the set of coupled channel equations and boundary conditions that determine the full 
S-matrix which is the focus of our study 



(T(r)-£I + y(r,A))T(r,A)=0, 



(2.12) 



T(r,A) ^ ^ ( h-{kr)-ifi+{kr)S{K\) 



(2.13) 
(2.14) 



The matrices (kr) are diagonal matrices with the incoming and outgoing spherical 
Riccati-Hankel functions on the diagonal 



h (kr)) = Si,h^ (kr), (h-^ (kr)) = S^h^ (kr) 



(2.15) 



The multi-channel S-matrix S is the object of interest that we will use to find the resonant 
states. We can extract the S-matrix from the numerical solution of (2.12) with its boundary 
conditions by evaluating the Wronskian in a point r outside the range of the potential 



S{k,A) = W Y{r,A),h+{kr) 



Yir,A),h-{kr) 



(2.16) 



where W denotes the Wronskian of two matrices defined as 



(2.17) 



with A, Be C"^". 
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3 Resonances and regularization 

The aim of this article is to develop a robust method to trace the parameter dependence of 
resonances that appear in multichannel problems. We define a resonance or bound state 
to be the solution of equation (2.1) for an energy E G C for which the S-matrix has a pole. 
This definition is widely used and it is motivated, amongst others, in [ , ]. 

It is natural to locate the resonances numerically by searching for the zeros of S{k,\)^^. 
However, it is not the best numerical strategy to search for these zeros directly. A key 
property of the S-matrix is the symmetry S{k*,\)=S{k,\)^^, where k* denotes the complex 
conjugate of k. This makes S{k,A) unitary for A: G R. As a result a pole of the S-matrix at k 
is always accompanied by a zero of S{k,A) at k*. This makes the search for the zeros of 
S{k,A)~^ a numerically hard problem, particularly when k<^l. In that case S{k,A) has 
zeros and poles close to each other and it does not satisfy the smoothness conditions 
required to ensure rapid convergence if, for instance, Newton iterations are used to find 
the zeros. 

For single-channel radial problems, as in section 2.1, it was proposed in [ ] to use a 
regularized function rather than the S-matrix itself to converge and to track down bound 
states and resonances. This function was defined as, per /-channel, 

with Si as in (2.4). It was shown that for k^O this function has a zero iff S/ {k,A) has a pole. 
In addition this function is bounded as long as drii{kr)V{r,A)ipi{k,r) ^^0, which does 
not hold only for very specific potentials. 

In this paper we extend this result to the more general class of coupled channel 
problems described in the previous section. We introduce, in extension of (3.1), the 
(^max + l) by (/max + 1) Square matrix function 

FiKA) = IC{S{KA)-I)-\ (3.2) 

where AT is a diagonal matrix with elements {)C)ii = k^^'^^ on the diagonal and where we 
extract the S-matrix using (2.16) from the numerical solution of (2.12), (2.13) and (2.14). 
The fimction used for continuation then is 



det(F(fc,A))= n^^'^^ / det{S{k,A)-I). (3.3) 




In appendix B we show that det(F(A:,A)) has a zero if and only if det(S(A:,A)^^) has 
a zero. This means that we can look for the bound states and resonances via (3.3) and 
investigate the application of continuation to equation det(F(fc,A)) = . In addition, the 
coalescing of zeros and poles that occurs in the S-matrix is now absent. Unfortunately, we 
cannot derive a bound for det(F(A:,A)) itself as was possible for (3.1). See appendix B for 
details. 
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4 Numerical Continuation methods 

In the study of dynamical systems one is often interested in the dependence and evolution 
of solutions of nonlinear problems in terms of some system parameter. Generally these 
solution sets can only be approximated numerically and finding those requires efficient 
and robust techniques. More specifically, numerical continuation techniques are con- 
cerned with approximating the solution set {x|F(x) = 0} of an underdetermined system 
of nonlinear equations: 



that depends on some system parameter A. Such equations arise for instance in the study 
of dynamical systems or parametrized systems of ODEs but virtually every problem that 
can be written in the above formulation and that satisfies certain smoothness conditions 
can be subjected to those methods. 

Due to the computationally intensive nature of many of these problems, efficiency is 
one of the main concerns in the numerical study of such dynamical systems. The number 
of evaluations of the function F(x) should therefore be kept to a minimum. Additionally, 
solution sets often exhibit complex geometries with intersecting and bifurcating solution 
branches. Generally, the study of bifurcating branches involves rigorous stability analysis 
of the underlying solutions and is a complicated subject on its own. In this work the 
only type of bifurcations that can occur, are so called "simple" bifurcation points. These 
manifest themselves as two intersecting solution branches and are characterized by the 
dimension of the null-space of the Jacobian Fx{xt) being 2 in a point Xt, which is called a 
"branch point". In many problems, however, the bifurcations are much more involved and 
a thorough treatise on bifurcations and stability of solutions can be found in [ - ]. 

In the last few decades various techniques for numerical continuation have been 
developed. Many of those fall in the category of so called "predictor-corrector" methods: 
starting from a known solution point Xj (such that F{xi) = 0), the next one is found by first 
taking a guess based on the previous point and correcting that guess in a second step. 
As such, starting from one known initial point xq in the solution set an approximation 
{xi\i = 0,...,n with F(x,) =0} of that solution set is constructed by finding subsequent 
points on the curve. 

Many choices can be made for the implementation of both steps. A good overview can 
be found in [ - ]. The popular "pseudo arclength" method [ ] takes a predictor-corrector 
approach. The predictor step uses the tangent t,- to the solution curve at to construct a 
prediction xf^-^ for the next point Xj+i. The corrector step consists of Newton iterations 
on the system of equations F = augmented with an additional equation that constrains 
the corrector step to a plane orthogonal to the tangent t,. The Newton iterations then 
converge to the intersection of that plane and the solution curve, guaranteeing a unique 
next point Xj+i in the continuation [10]. 




(4.1) 
(4.2) 
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The pseudo arclength algorithm has been implemented in the AUTO [11, 12] library 
which we have used in this work. Other well-known implementations of numerical 
continuation algorithms include LOCA [ ] which is a part of the Trilinos framework [ ' 4], 
Matcont [15] (a Matlab implementation) and Multifario [16] which allows multi-para- 
meter continuation. 

5 Numerical Results 

The results of this section have been obtained with a renormalized Numerov solver [17,18] 
for the calculation of the coupled channel solutions of equation (2.12). However, these 
results are independent of the particular numerical method used to solve the equations as 
long as the method is able to extract the S-matrix e.g. through equation (2.16). Indeed, the 
regularized function det(F(A:,A)) is the only required input for the continuation routines. 

To demonstrate our approach we have applied numerical continuation to a number 
of model problems involving well-behaved short-range potentials. The numerical im- 
plementation of the continuation algorithm itself is provided by the AUTO-07p Fortran 
library [ ]. The coupled channel renormalized Numerov solver and the code to obtain 
the regularized function det(F(A;,A)) are used as driver routines for AUTO and we have 
implemented them in C++. 

5.1 Comparison with analytically solvable models 

As a first testing example for the study of numerical code and regularization we have 
compared our results thoroughly with a model problem from [ ]. There, among other 
examples, a short-range two-channel coupling model for the low-energy electron-COa 
scattering is studied. An analytically solvable expression with square-well potentials is 
presented together with resulting bound/ resonant pole trajectories. 

The application of numerical continuation on both the analytical expressions and 
the numerical solutions of this model matched the results previously published by the 
authors. 

5.2 Gauss potential with s-wave and p-wave coupling 

In this section we look at a system of two coupled Gaussian potential wells with angular 
momenta /i =0 (referred to as s-wave in physics) and /a = 1 (referred to as p-wave). The 
mass is chosen to be p = 1. The two channel potentials and the coupling potential are set 
to: 

Vij{r,\) = -k,je-''\ 

where z,;G {1,2} and Vij{r,X) = Vji{r,X). 

Let us first consider each of the two channels separately, i.e. uncoupled (A12 = 0). 
Setting All, the strength of the first-channel potential, to 7 gives a system with one bound 
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state shown in table 1(a). Similarly, bound states in the second channel with A22 = 20 
are shown in table 1(b). These results were calculated through Newton iterations on the 



(a) First channel. An =7 (b) Second channel, A22 =20 



i 


ki (Newton) 


ki (Matslise) 


i 


ki (Newton) 


ki (Matslise) 





2.185562e+00i 


2.185562e+00i 





3.617543e+00i 


3.617478e+00i 








1 


8.938842e-01i 


8.938842e-01i 



Table 1: Bound states of uncoupled channels in the s- and p-wave example, obtained with (1) Newton iterations 
on det(f (fc,A)) (Newton) and (2) with Matslise. 

function det(F(A;,A)). They compare, see table 1, favorably with results obtained with 
Matslise [20]. 



i 


ki (Newton) 


ki (Matscs) 


ki (fin. DIFF.) 





3 


623677e+00i 


3 


623551e+00i 


3 


.623677e+00i 


1 


2 


178012e+00i 


2 


178012e+00i 


2 


. 195645e+00i 


2 


9 


035406e-01i 


9 


035275e-01i 


9 


.039594e-01i 



Table 2: Coupled channel bound states in the s- and p-wave example, with Newton iteration (Newton), coupled 
channel extension of CPM{P,N} (Matscs) and finite differences (fin. diff.) 

The coupled system on the other hand has bound states at positions deviating from 
those in table 1 as the coupling strength increases. In this case, with A12 = 0.5, our results, 
see table 2, are in close agreement with those obtained through the recently extended 
CPM{P,N} method for coupled channel systems [21,22] implemented in MATSCS [ ]. 
Reference values obtained with eigenvalue calculations based on finite difference dis- 
cretization, although not very accurate, are also provided. 

It is now interesting to see how these bound states evolve in terms of variations in A22. 
Applying numerical continuation with starting values fc, (Newton) taken from table 2 
gives results shown in figure 1. The continuation was performed both in the direction of 
the origin and away to show as much of the resonant trajectories as possible. 

Details of the underlying numerical procedures are as follows. The integration interval 
of the coupled channel Numerov method spans from to 4.6 and contains 4096 points. 
The tolerance for Newton iterations used by AUTO is of the order 10^^ whereas the length 
of the continuation steps might vary adaptively between lO^'^ and 10^^. 

Each of the three continuation branches originates from a bound state and follows 
a path towards the resonant region (i.e. the lower half of the complex /c-plane) as the 
potential strength A22 is decreased. Near the origin of the complex plane the second 
channel potential becomes too weak to hold the bound state and reaches a threshold 
value where the state transforms into a resonance. Looking at the continuation curves 
in figure 1 we can therefore distinguish three regions: the bound state regime on the 
positive imaginary axis, the transition into a resonance around the origin and the resonant 
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regime below the real axis. Their properties are summarized as follows. (1) The bound 
states are situated on the positive imaginary axis and correspond to real, negative (rel. to 
potential asymptote) energies. In the neighborhood of the bound states of the first channel, 
in this case k ^ 2.2z, "avoided crossings" can be seen where a relatively big decrease in 
A22 has a minor influence on kj. However, as soon as the lower state approaches, 
the higher one gets pushed away from this position and makes place for (2) As 
A22 passes certain threshold values, the potential becomes too weak to hold a bound 
state ki which transitions into a resonance. For s-waves these transitions happen on the 
negative imaginary axis not far from the origin, whereas for higher angular momenta this 
branching occurs exactly at 0. In the case of s/ p-wave coupling, both effects influence each 
other depending on the coupling strength and bound states of the separate channels. The 
threshold values for A22 with the corresponding branch points for the three states in our 
example are summarized in table 3. (3) In the resonant regime (Im(A;) < 0) the states gain a 
real part which corresponds to a positive scattering energy. The negative imaginary part 
of these states is related to the time delay in the time-dependent picture. These branches 
are symmetric but only those in the fourth quadrant have physical significance. On the 
negative imaginary axis one also finds so called virtual states. 



i 


Reiki) 


lm{ki) 




1 


-1.227997e-12 


-7.999927e-04 


6.091215e+00 


2 


6.280509e-13 


-1.589331e-02 


1.742094e+01 



Table 3: Branch points at which bound states transition into a resonance for the Gauss s/p-wave example. The 
real part of k is negligible whereas the imaginary part indicates slight off-zero branching discussed in 5.2. 

The proposed regularization (3.3) plays a significant role near and at threshold values. 
Indeed, we have proven in appendix B the absence of additional, spurious zeros of 
det(F(A:,A)). All A22-dependent poles (which would otherwise coalesce with zeros in the 
origin at threshold values) are eliminated as well which makes continuation through 
the origin possible. We can not, however, guarantee smoothness of det(F(A;,A)) over the 
whole domain as demonstrated by the two stationary poles in figure 2, where we show a 
comparison of det(S(A:,A)) and det(F(A:,A)). 

5.3 Gauss potential with p-wave and rf-wave coupling 

The second model problem also uses Gauss potentials in both channels as well as the 
coupling and ji = l. However, angular momenta for channel 1 and 2 were chosen to be 
li = l and /2 =2. This choice is motivated by the previously mentioned distinction between 
s-waves and waves with higher angular momenta in 5.2: In the / = case the bifurcation 
does not take place in the origin of the /c-plane but on the imaginary axis below the origin, 
whereas for / > it does branch off in the origin [ ]. By taking both angular momenta 
higher than zero we enforce branching of the coupled channel system in the origin and 
therefore an essentially different geometry for testing the regularized formulation than in 
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the previous example. 

Setting All = 10, A22 = 30 and A12 = 0.3 gives a system with 3 bound states. The 
approximations, analogous to the first example, are given in table 4. 



i 


ki (Newton) 


ki (MatSCS) 


ki (fin. DIFF.) 





3 


.796532e+00i 


3 


796472e+00i 


3 


.796532e+00i 


1 


1 


.600083e+00i 


1 


600083e+00i 


1 


.600127e+00i 


2 


6 


.599123e-01i 


6 


599112e-01i 


6 


.603354e-01i 



Table 4: Bound states of the coupled channel p- and d-wave example. We compare values obtained by Newton 
iterations (Newton) on det(F(fc,A)), (Matscs) and finite differences (fin. diff.) 

Again, these values were used as starting points for numerical continuation. The 
resulting curves are shown in figure 4. The details of the coupled channel renormalized 
Numerov solver and AUTO continuation parameters were set to the same values as in the 
first example. 

Also in this case the proposed regularization provides a viable way of removing poles 
in the origin at thresholds, hereby allowing continuation of the relevant zeros through the 
origin. This is shown in figure 3. This figure also shows the sharpness of nearby poles 
and zeros of det(S(A:,A)) compared to much smoother behavior of det(F(A:,A)) which 
positively influences the convergence of Newton iterations in the continuation process. 

6 Discussion and Conclusions 

In this article we have applied numerical continuation to track resonant states in a coupled 
channel Schodinger equation. In our approach, we track the zeros of a function that 
can be extracted from the numerical solutions of the coupled channel equation. This 
function, related to the S-matrix by a transformation, has a zero wherever the S-matrix has 
a pole and is amenable to numerical continuation. The analysis has also revealed that the 
function can have accidental poles which may slow down the convergence of the Newton 
iterations. However, for our model problems, we have not seen any difficulties caused by 
this. 

We have tested this approach and solved a range of two channel model problems with 
robust results. The method automatically tracked all the resonant states and detected the 
bifurcation points. It also successfully traversed the threshold which indicates a successful 
regularization. 

Our motivation to develop this method originates from chemical reactions that are 
modeled by the Born-Oppenheimer approximation. It is based on slow-fast separation 
between the motion of nuclei and electrons that is justified because of their mass difference. 
In this approximation the modeling is a two step process. In the first step the scattering 
of the electrons in the field of the nuclei is solved with the positions of the nuclei as 
parameter in this scattering problem. It is especially important to identify the resonances 
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in the scattering problem. In the second step the molecular dynamics of the nuclei, that 
move on the potential surface given by the resonances determined in the first step, is 
solved. 

This paper is focused on the first step where the resonances need to be identified in a 
family of scattering problems each with a different parameter, which is the position of the 
nuclei in these problems. 

In our future work we aim to create a robust and reliable tool that scales to realistic 
systems. 
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A Spherical functions 

Various solutions of the free radial wave equation are intensely used throughout this 
paper as they can describe the asymptotic behavior of the wave functions that solve the 
(non-free) Schrodinger equation. In an attempt to avoid ambiguities we summarize the 
used conventions in this appendix. For a far more thorough treatment we suggest [24] 
and [ ], where similar conventions are used. 

Spherical Riccati-Bessel functions are defined as solutions of the free radial Schrodinger 
equation 





(A.1) 



The different, commonly used particular solutions of this equation are: {x = kr) 



• ji{x) = xjj{x) = iyx^//+i/2(x), regular solutions, spherical Riccati-Bessel functions. 




functions. E.g.: no(^) = — cos(x), ni(x) = — ^cos(x) — sin(x), ... 



• hfix) = —n/(x)±zy/(x) Incoming (—) and outgoing (-|-) spherical Riccati-Hankel 




B Regularization lemma 



This appendix technically describes the proof that the application of the regularizing pro- 
cedure in section 3 does not introduce any zeros of the continuation function det(F(A:,A)) 
that can not be attributed to bound or resonant states. The proof involves the use of 
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some well-known definitions and results from quantum mechanical scattering and gives 
a thorough explanation of the procedure yet is not necessary for the understanding and 
reproduction of the main results of our contribution. The conventions and definitions 
used are mostly those foimd in [ ] and, on occasion in [5,25]. Section B.l summarizes 
briefly a necessary minimal subset of definitions for the proof whereas section B.2 proves 
why the regularizing procedure does not introduce false resonances or bound states. 

For clarity we omit the potential parameter A in our notations as it does not play a role 
in the derivations; unless it is really necessary or clarifying. Following [4] we do however, 
indicate the explicit dependence of the wave functions onk,e.g.Y{k,r). 



B.l The Jost matrix and the S-matrix 



The "regular solution" 0{k,r) of equation (2.12) is defined to behave exactly as }{kr) at 
r = 0. This fixes both ^ and its normalization uniquely. At infinity the regular solution 
must tend to a linear combination of two other linearly independent solutions such as the 
free solutions h^, hence allowing us to write 



0iKr) ^^[h- {kr)f{k) - h+ {kr)f{ky 



(B.l) 



which defines the "Jost matrix" f{k) as a diagonal matrix that bundles so called "Jost 
functions" fi{k) 

(m \ 



\ 



(B.2) 



(it is known that fi{k)* =fi{—k) and hence f{k)* =f{—k)) Since we have already written 
the asymptotic form of the "physical solution" Y{k,r) (which is calculated by the numerical 
solver) in terms of the S-matrix as in 



Y{Kr) 



''^%llh-{kr)-h+{kr)S{k)), 



it follows that the S-matrix can be expressed as the ratio of the Jost matrices 

S{k)=f{-k)f{k)-' 
and that and Y are related through the well-known expression^ 

Y{Kr) = 0{Kr)f{k)-\ 



(B.3) 



(B.4) 



(B.5) 



Please note the different approaches studied in various reference works on the subject. Some authors such 
as [ ] prefer the simpler approach of expression (B.5), whereas others like [ ] explicitly write certain factors 
without incorporating them, for instance, in the Jost function. In the single channel variant of (B.5) for instance 

this would result irnpi{k,r)= (^2i+i)n '^fi{k)^ • Unless explicitly stated otherwise, we choose to follow the former 
convention. 
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The previous equations illustrate a widely accepted result of bound and resonant states 
being represented, equivalently, by poles of the S-matrix and zeros of the Jost function. As 
we can not extract the Jost matrix directly from the numerical solution of (2.12) but only a 
function that is proportional the Jost matrix, we have to use the poles of the S-matrix as 
the characterization of bound and resonant states. 

The physical solution Y can also be written as the solution of the coupled channel 
partial wave version of the Lippmann-Sch winger equation 

/"OO 

Y{r)=J{kr)+ dr'G{r,r'-k)U{r')Y{r'), (B.6) 

JO 

where G{r,r';k) denotes the diagonal matrix of one-charmel Green's functions G^ {r,r';k) = 

^;;(A;r<)^^(A;r>), with r< =min(r,r') and r> =max(r,r') and where U = 2}iV. If we now 
combine (B.6) with the asymptotic behavior of (2.14) we find 

S(k)-I=^ dr'j(kr')Uir')Y(r'). (B.7) 
k Jo 

Inside the integral we replace Y with the regular solution through the identity (B.5) 
which leads to the equation 

S{k)-I = =^n"^dr'jXkr')U{r')0{r'))f{kr\ (B.8) 



B.2 Bounds on the regularization 

Equation (B.8) states that the (/max + 1) x ('max + 1) matrix S{k) — I is proportional to a 
matrix of integrals multiplied by the inverse of the Jost matrix. To guarantee that the 
regularizing procedure from section 3 defined by eq. (3.3) as 

('max \ / 

Ylk^'+'U det{S{KA)-I), (B.9) 

does not introduce false bound states or resonances we must prove that zeros of the Jost 
matrix f{k) are the only contributors to the poles of expression (B.8) and as such, to the 
zeros of det(F(fc,A)). This requires the matrix of integrals in equation (B.8) to be bounded 
which is guaranteed by the following lemma. 

An important remark has to be made ragarding the regularity constraints of the 
potential matrix V. The proof of lemma B.l is a generalization of the one charmel case 
where coupling terms in the potential do not occur. There, the boimdary conditions that 
define the regular solution fi and the iterative procedure of solving the one channel 
analogue of integral equation (B.6) for cpi indeed would give rise to proper bounds. 
Contrary to the one channel case however, the off-diagonal terms in V that couple channels 
with different angular momenta result in serious singularities in the origin r = 0, hereby 
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destroing the possibility for a bound without very restrictive assumptions on the potentials 
that would cancel out the singular behavior in the origin [ , ]. These restrictions can be 
lifted to the much more liberal requirement for the off-diagonal terms drr\Vi^j{r) \ < oo 
by rewriting the integral equation (B.6) with additional terms. The iterative procedure to 
solve the integral equation then gives the same regularity properties of as in the one 
channel case. 

To retain clarity we demonstrate the bound using the more restrictive assumptions for 
the off -diagonal potential terms, keeping in mind that a similar bound can be found in 
the more liberal case of, for instance, the potentials studied in the numerical examples in 
section 5. 

Lemma B.l. Let be the regular solution of eq. (2.12) where U = 2piV the potential matrix 
satisfies the (restrictive) short range conditions 

f/~drr'.+'.+^|T^.,(r)|<co V.,; 

\/~drr-''+'.+i|y,(r)|<co Vz,;, ^''•''^^ 

and let jikr) he the diagonal matrix of spherical Riccati-Bessel functions with corresponding 
then the integral 

1= / dri(kr)U(r)0(r), (B.ll) 

JO 

zs bounded. 

Proof. Because of the particular boundary conditions for r — > and the restrictions on the 
potential matrix V the regular solution fits the matrix integral equation (20.19) in [4] 

0{r) = i{kr)+ f dr'g(ry)U{r')0{r'), (B.12) 

JO 

where g{r,r') is a diagonal matrix composed of Green's functions 

gi, (r,r') = (1/fc) (/,, {kr)^ {kr') - n,, {kr)Ji^ (/cr')) , (B.13) 

on the diagonal. 

This regular solution matrix can be written as a series 0{r) =J^„ ^'"^ (r) with the first 
term 0^^^ =j{kr) and the following recurrence relation between 0^"^ and 0^"^^^ 

0("){r)= f dr'g{ry)U{r')0^''-^\r'). (B.14) 

JO 

The integral (B.ll) can then be written as 

oo ^oo 

X=Y dr]{kr)U{r)0^"\r) (B.15) 
= / drj{kr)U{r)j{kr)+ / drf{kr)U{r) / dr,g{rj,)U{ri)j{kri) (B.16) 

Jo JO JO 

H h / drj{kr)U{r) / drig{r,ri)U{ri) / dr2--- / drngir„^i,rn)U{rn)jikr„) + - ■ ■ 

Jo JO JO JO 



(B.17) 
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Without loss of generality, let us look for example at a 2 x 2 coupled channel problem 
with /i in the first channel and I2 in the second channel. The integral above is then 





/ dr 




Jo 






+ 


/ dr 




10 






X 


fdr, 




/o 













) 



U,2{r)\ (kikr) ^ 
U2i{r) U22ir)J\ kikr)^ 

Un{r) Ui2{r 
2i{r) U22{r 

Uuin) U^2{ri)\ (Jkikri) 
U2i{ri) U22{ri)J V 







(B.I8) 

(B.19) 
(B.20) 



When we look at a particular element of this two by two example we see the contribu- 
tions 



^11= drJi^{kr)Uii{r)Ji^{kr)+ drJi^{kr)Uii{r) dr^gi^{r ,r^)Un{r^)k{kr^) 

1/ "-'0 "-'0 

/•oo i-r 

+ dr]i^{kr)Ui2{r) j^dr^gi^{r,r^)U2i{r^)k{kr^) + .... (B.21) 



We note that Riccati-Bessel functions that appear in the integral both have I = h. In a 
similar way the matrix element lu will have the ji^ in the beginning of the integral and //^ 
at the end. 

To show that this integral Xn is finite we make use of the bounds fond in [4,5,25]. 
There exist constants Ci and C2 such that 



i;/(^r)|<Cieli-«l'1- 



\k\r 



^ + \k\r 



|g/(r,/)|<C2eli-«l('-^')^ 





\k\r 


1+1 


k\ 


\r 



/+1 



l + \k\r' 



(B.22) 
(B.23) 
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After substitution we find 

2/1+2 



+='r*^""'*(T^)""l""(')l=^'""'""(l^)"" 



X 



= A;^'^+^Cii, for some constant Ci i . 

This series can be summed up because the potentials Vn, V12, V21 and V22 satisfy the 
regularity condition (B.IO). In a similar way we can find bounds for the other elements of 
1. 

h\ a more general case we have the bound 



/ 12^11 1 IX12I . 

\^2l\ \^22\ I < 

V '■■ 



This indeed means that the only way expression (B.9) can have a zero is through a zero of 
the Jost function, which ensures that no additional false "resonances" or "bound states" 
are introduced by using zeros of this function to locate them. Naturally every zero of 
f{k,\) will result in a zero of det(F(A:,A)) which means we are able to track all resonances 
and bound states through the use of (B.9). This proves the equivalence of boimd and 
resonant states and the zeros of det(F(A;, A) ) . □ 
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(c) Continuation curves in full space with projections along the sides. 



Figure 1; Continuation curves of the Gauss s/p-wave coupling example. Branches of different colors originate 
from different starting points (drawn as black crosses) in table 2: i = (red), i = l (blue), i = 2 (green) 
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Figure 2: Contour plots of |det(S)| (upper row) 
and |det(f)| (lower row) in the Gauss s/p-wave 
example (logarithmic scaling was applied to "flat- 
ten" the images). Ai2 = 0.5, An =7 and A22 = 17.5 
(left) and 17.4 (right). As A22 is decreased, poles 
of the S-matrix move from the bound state and 
virtual state region through the origin of the com- 
plex fc-plane into the resonant regime {i = 2 state 
becomes resonant). Arrows indicate poles (P) and 
zeros (Z). 




Figure 3: Contour plots of |det(S)| (upper row) 
and |det(F)| (lower row) in the Gauss p/d-wave 
example (logarithmic scaling was applied to "flat- 
ten" the images). Ai2=0.3, An =10 and A22 = 13-5 
(left) and 13.4 (right). As A22 is decreased, poles 
of the S-matrix move from the bound state and 
virtual state region through the origin of the com- 
plex plane into the resonant regime (f = l state 
becomes resonant). Arrows indicate poles (P) and 
zeros (Z). 



22 





lm(k) 

Re(k) 

(c) Continuation curves in full space with projections along the sides. 

Figure 4: Continuation curves of the Gauss p/d-wave coupling example. Branches of different colors originate 
from different starting points (drawn as black crosses) in table 4: i = (red), i = l (blue), i = 2 (green) 



